home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / DBio / DSequence.cpp < prev    next >
Encoding:
Text File  |  1996-07-05  |  39.1 KB  |  1,685 lines  |  [TEXT/R*ch]
  1. // DSequence.cp 
  2. // d.g.gilbert
  3.  
  4.  
  5.  
  6. #include "DSequence.h"
  7. #include "DSeqFile.h"
  8. #include "DSeqList.h"
  9. #include "DREnzyme.h"
  10. #include <ncbi.h>
  11. #include <dgg.h>
  12. #include <DUtil.h>
  13. #include <DFile.h>
  14. #include "ureadseq.h"
  15.  
  16.  
  17. char* DSequence::kConsensus = "Consensus";
  18. char  DSequence::indelHard = '-';
  19. char  DSequence::indelSoft = '~';
  20. char  DSequence::indelEdge = '.';
  21.  
  22.  
  23. Local char* gAminos        =    "ABCDEFGHIKLMNPQRSTVWXYZ*_.-~?";
  24. Local char* gIubbase    = "ACGTUMRWSYKVHDBXN_.-*~?";     //did include ".", for what?
  25. Local char* gPrimenuc    = "ACGTU_.-*~?";
  26. Local char* protonly    = "EFIPQZ";
  27. Local char* stdsymbols= "_.-*?";
  28. Local char* allsymbols= "_.-*?!=/**/[]()!&#$%^&=+;:/|`~'\"\\";
  29. Local char* seqsymbols= allsymbols;
  30.  
  31.  
  32.  
  33.  
  34.  
  35.  
  36. // DSequence statics -----------------
  37.  
  38.  
  39. long DSequence::NucleicBits(char nuc)
  40. {
  41.     //char nuc= toupper(nuc);
  42.     switch (nuc) {
  43.         case 'a': case 'A': return kMaskA;  
  44.         case 'c': case 'C': return kMaskC;   
  45.         case 'g': case 'G': return kMaskG;   
  46.         case 't': case 'T': return kMaskT;   
  47.         case 'u': case 'U': return kMaskU;  // ==  ( kMaskT | kMaskExtra); 
  48.         case 'y': case 'Y': return (kMaskC | kMaskT);  
  49.         case 'k': case 'K': return (kMaskG | kMaskT);  
  50.         case 's': case 'S': return (kMaskG | kMaskC);  
  51.         case 'r': case 'R': return (kMaskG | kMaskA);  
  52.         case 'm': case 'M': return (kMaskA | kMaskC);  
  53.         case 'w': case 'W': return (kMaskA | kMaskT);  
  54.         case 'h': case 'H': return (kMaskA | kMaskC | kMaskT);  
  55.         case 'b': case 'B': return (kMaskG | kMaskC | kMaskT);  
  56.         case 'v': case 'V': return (kMaskG | kMaskC | kMaskA);  
  57.         case 'd': case 'D': return (kMaskG | kMaskT | kMaskA);  
  58.         case ' ': case '_': return 0;  //? spacers
  59.         case  0 : return 0;  
  60.         case 'n': case 'N': return kMaskNucs; 
  61.         default : 
  62.             {
  63.             if (nuc == indelHard || nuc == indelSoft)  return kMaskIndel;
  64.             else if (nuc == indelEdge) return 0; //return kMaskNucs; ???
  65.             else return 0; //kMaskNucs;  //?? match all or match none ??
  66.             }
  67.         }
  68. }
  69.  
  70.  
  71. Boolean DSequence::IsNucleicMatch( long xNucBits, long yNucBits)
  72. {
  73.         // note: Indel/spaces match nothing, including other indel/space ...
  74.     return 0 != ( (kMaskNucs & xNucBits) & (kMaskNucs & yNucBits));
  75. }
  76.  
  77.  
  78. char DSequence::NucleicConsensus( long xNucBits, long yNucBits)
  79. { 
  80.     char nuc; 
  81.     Boolean isrna, indel;
  82.  
  83.     if (xNucBits==0 || yNucBits==0) 
  84.         return ' ';  //?? either a space yields space ?
  85.     else {
  86.         indel= (kMaskIndel == xNucBits || kMaskIndel == yNucBits);
  87.         isrna= (kMaskU == xNucBits || kMaskU == yNucBits);
  88.         switch ( (kMaskNucs & xNucBits) | (kMaskNucs & yNucBits)) {
  89.             case kMaskA: nuc= 'A'; break;
  90.             case kMaskC: nuc= 'C'; break;
  91.             case kMaskG: nuc= 'G'; break;
  92.             case kMaskT: 
  93.                 if (isrna) nuc= 'U'; else nuc= 'T';
  94.                 break;
  95.                     //? do other ambig codes ?
  96.             default:    nuc= '.'; break; //? use this as Consensus other ?
  97.             }
  98.         if (indel) nuc= tolower(nuc); //?
  99.         }
  100.     return nuc;
  101. }
  102.  
  103.  
  104.  
  105. void DSequence::NucleicComplement( Boolean isrna, char* fromSeq, char* toSeq, long len)
  106. {
  107. // note: fromSeq may == toSeq 
  108.     register char b, c;
  109.     Boolean islow;
  110.      while (len--) {
  111.         b= *fromSeq++;
  112.         islow= islower(b);
  113.         if (islow) b= toupper(b);
  114.         switch (b) {
  115.             case 'C': c= 'G'; break;
  116.             case 'G': c= 'C'; break;
  117.             case 'A': if (isrna) c= 'U'; else c= 'T'; break;
  118.             case 'T':
  119.             case 'U': c= 'A'; break;
  120.             case 'R': c= 'Y'; break;
  121.             case 'Y': c= 'R'; break;
  122.             case 'M': c= 'K'; break;
  123.             case 'K': c= 'M'; break;
  124.             case 'S': c= 'S'; break;  // was mistaken 'w'
  125.             case 'W': c= 'W'; break;  // was mistaken 's'
  126.             case 'H': c= 'D'; break;
  127.             case 'D': c= 'H'; break;
  128.             case 'B': c= 'V'; break;
  129.             case 'V': c= 'B'; break;
  130.             default : c= b;   break;
  131.             }
  132.         if (islow) c= tolower(c);
  133.         *toSeq++= c;
  134.         }
  135. }
  136.  
  137.  
  138.  
  139. const char* DSequence::Amino123( char amino1)
  140. {
  141.  switch (toupper( amino1)) {
  142.      case 'A': return "Ala";
  143.      case 'M': return "Met";
  144.      case 'R': return "Arg";
  145.      case 'F': return "Phe";
  146.      case 'N': return "Asn";
  147.      case 'P': return "Pro";
  148.      case 'D': return "Asp";
  149.      case 'S': return "Ser";
  150.      case 'C': return "Cys";
  151.      case 'T': return "Thr";
  152.      case 'Q': return "Gln";
  153.      case 'W': return "Trp";
  154.      case 'E': return "Glu";
  155.      case 'Y': return "Tyr";
  156.      case 'G': return "Gly";
  157.      case 'V': return "Val";
  158.      case 'H': return "His";
  159.      case 'I': return "Ile";
  160.      case 'B': return "Asx";
  161.      case 'L': return "Leu";
  162.      case 'Z': return "Glx";
  163.      case 'K': return "Lys";
  164.      case 'X': return "Xaa";
  165.      case '*': return "End";
  166.      case ' ': return "   ";
  167.      default : return "???"; //??
  168.      }
  169. }
  170.  
  171.  
  172. char DSequence::Amino321( char* amino)
  173. {
  174.     char    aa[4];
  175.     aa[0]= toupper(amino[0]);
  176.     aa[1]= tolower(amino[1]);
  177.     aa[2]= tolower(amino[2]);
  178.     aa[3]= 0;
  179.     
  180.              if (!StrCmp(aa,"Ala")) return 'A';
  181.     else if (!StrCmp(aa,"Met")) return 'M';
  182.     else if (!StrCmp(aa,"Arg")) return 'R';         
  183.     else if (!StrCmp(aa,"Phe")) return 'F';
  184.     else if (!StrCmp(aa,"Asn")) return 'N';       
  185.     else if (!StrCmp(aa,"Pro")) return 'P';
  186.     else if (!StrCmp(aa,"Asp")) return 'D';    
  187.     else if (!StrCmp(aa,"Ser")) return 'S';
  188.     else if (!StrCmp(aa,"Cys")) return 'C';         
  189.     else if (!StrCmp(aa,"Thr")) return 'T';
  190.     else if (!StrCmp(aa,"Gln")) return 'Q';       
  191.     else if (!StrCmp(aa,"Trp")) return 'W';
  192.     else if (!StrCmp(aa,"Glu")) return 'E';    
  193.     else if (!StrCmp(aa,"Tyr")) return 'Y';
  194.     else if (!StrCmp(aa,"Gly")) return 'G';          
  195.     else if (!StrCmp(aa,"Val")) return 'V';
  196.     else if (!StrCmp(aa,"His")) return 'H';
  197.     else if (!StrCmp(aa,"Ile")) return 'I';       
  198.     else if (!StrCmp(aa,"Asx")) return 'B';
  199.     else if (!StrCmp(aa,"Leu")) return 'L';          
  200.     else if (!StrCmp(aa,"Glx")) return 'Z';
  201.     else if (!StrCmp(aa,"Lys")) return 'K';           
  202.     else if (!StrCmp(aa,"End")) return '*';           
  203.     else if (!StrCmp(aa,"   ")) return ' '; //??           
  204.     else return 'X';
  205. }
  206.  
  207.  
  208. char* DSequence::Nucleic23Protein( char* nucs, long nbases)
  209. //!! Need to do this w/ NucBits so ambiguity codes work !!
  210. //!! The CodonTable.codon == acodon should be IsNucleicMatch( bits(codon), bits(acodon)) 
  211. {
  212.     long    i;  //borland 16bit don't like longs for iterators !?
  213.     char  *ap, *kp, *aminos;
  214.  
  215.     if (DCodons::NotAvailable()) return NULL; // (gCodonTable[0].amino == 0)
  216.  
  217.     aminos= StrDup( nucs);
  218.     for (i= 0, ap= aminos; i<nbases; i++) {
  219.         register char c= toupper( *ap);
  220.         if (c=='U') c= 'T';
  221.         *ap++= c;
  222.         }
  223.  
  224.     for (i= 0, ap= aminos, kp= aminos; i<nbases; i++) {
  225.         for (long j=0; j<64; j++)
  226.             if ( MemCmp(DCodons::codon(j), ap, 3) == 0) { //gCodonTable[j].codon
  227.                 *kp++= DCodons::amino(j); // gCodonTable[j].amino;
  228.                 break;
  229.                 }
  230.         }
  231.     *kp= '\0';
  232.     aminos= (char*) MemMore( aminos, kp - aminos + 1);
  233.     return aminos;
  234. } 
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244. // DSequence -----------------
  245.  
  246.  
  247. DSequence::DSequence()
  248. {
  249.     fBases= (char*) MemGet(1,true); 
  250.   fMasks= NULL; //(char*) MemGet(1,true); 
  251.   fMasksOk= false;
  252.     fInfo    = (char*) MemGet(1,true);     
  253.     fName = StrDup("Untitled");         //!? need unique name, eg. Noname-1,noname-2...
  254.     fKind = kOtherSeq;         //? user default
  255.     fLength        = 0;
  256.     fChecksum    = 0;
  257.     fIndex        = 0; //dummy
  258.     fModTime    = gUtil->time();
  259.     fChanged    = FALSE; //TRUE;
  260.     fIsMega        = FALSE;
  261.     fShowORF    = FALSE;
  262.     fShowRE        = FALSE;
  263.     fShowTrace = FALSE;
  264.  
  265.     fDeleted    = FALSE; 
  266.     fOpen            = FALSE;
  267.     fSearchRec.targ= NULL;
  268.     fOrigin    = 0;
  269.     fSelStart= 0;
  270.     fSelBases= 0;
  271. }
  272.  
  273.  
  274. DSequence::~DSequence()  
  275. {
  276.     MemFree(fBases);
  277.     MemFree(fMasks);
  278.     MemFree(fInfo);
  279.     MemFree(fName);
  280.     MemFree(fSearchRec.targ);
  281. }
  282.  
  283.  
  284. DObject* DSequence::Clone() // override 
  285. {    
  286.     DSequence* aSeq= (DSequence*) DObject::Clone();
  287.     aSeq->fBases= NULL; aSeq->SetBases( fBases, true);  
  288.     aSeq->fMasks= NULL; aSeq->SetMasks( fMasks, true); 
  289.     aSeq->fInfo=  StrDup( fInfo);
  290.     aSeq->fName=  StrDup( fName);
  291.     aSeq->fSearchRec.targ= StrDup( fSearchRec.targ);
  292.     return aSeq;
  293. }
  294.  
  295. void DSequence::CopyContents( DSequence* fromSeq)
  296. {
  297.     SetBases( fromSeq->fBases, true); 
  298.     SetMasks( fromSeq->fMasks, true); 
  299.     SetInfo( fromSeq->fInfo); 
  300.     SetName( fromSeq->fName); 
  301.     fKind            = fromSeq->fKind;          
  302.     fLength        = fromSeq->fLength;
  303.     fChecksum    = fromSeq->fChecksum;
  304.     fModTime    = fromSeq->fModTime;
  305.     fChanged    = fromSeq->fChanged;
  306.     fDeleted    = fromSeq->fDeleted; 
  307.     fOpen            = fromSeq->fOpen;
  308.     fIndex        = fromSeq->fIndex;
  309.     fSearchRec= fromSeq->fSearchRec;
  310.     fSearchRec.targ= StrDup( fromSeq->fSearchRec.targ);
  311. }
  312.  
  313. void DSequence::SetBases( char*& theBases, Boolean duplicate)
  314. {
  315.     long  newlen= (theBases) ? StrLen(theBases) : 0;
  316.     if (fBases) { 
  317.         if (fLength>0 && 
  318.             (fLength != newlen || MemCmp( theBases, fBases, newlen)!= 0))
  319.                 fChanged= TRUE;
  320.         MemFree(fBases); 
  321.         fBases= NULL; fLength= 0; 
  322.         }
  323.     if (theBases) {
  324.         if (duplicate) fBases = StrDup(theBases); 
  325.         else { fBases= theBases; theBases= NULL; }
  326.         fLength= newlen;  
  327.         }
  328.     fModTime= gUtil->time();
  329.     fMasksOk= false; // assume SetBases always called before SetMask...
  330. }
  331.  
  332.  
  333. void DSequence::UpdateBases( char* theBases, long theLength)
  334. {
  335.     fBases =  theBases;  
  336.     fLength= theLength;
  337.     fMasksOk= false; // assume this allways called before SetMask...
  338. }
  339.  
  340.  
  341. void DSequence::SetMasks( char*& theMasks, Boolean duplicate)
  342. {
  343.     if (fMasks) MemFree(fMasks);
  344.     if (theMasks) {
  345.         if (duplicate) fMasks = StrDup(theMasks);  
  346.         else { fMasks= theMasks; theMasks= NULL; }
  347.         FixMasks();
  348.         }
  349.     else {
  350.         fMasks= NULL;
  351.         fMasksOk= false;
  352.         }
  353. }
  354.  
  355. void DSequence::FixMasks()
  356. {
  357.     ulong i, mlen= 0;
  358.     if (!fMasks) { 
  359.         mlen= 1+fLength;
  360.         fMasks= (char*) MemGet(mlen,true);  
  361.         for (i=0; i<mlen; i++) fMasks[i] |= 0x80; // set high bit so Strxxx sees it as non-null
  362.         fMasks[mlen-1]= 0;
  363.         }
  364.     else { 
  365.         mlen= 1+StrLen(fMasks);
  366.         if (fLength+1 != mlen) {
  367.             fMasks= (char*) Nlm_MemExtend( fMasks, fLength+1, mlen);
  368.             mlen= fLength+1; 
  369.             for (i=0; i<mlen; i++) fMasks[i] |= 0x80; 
  370.             fMasks[mlen-1]= 0;
  371.             } 
  372.         }
  373.     fMasksOk= true;
  374. }
  375.  
  376. short DSequence::MaskAt(long baseindex, short masklevel)
  377. {
  378.     short b = kMaskEmpty;
  379.     if (fMasks && fMasksOk && baseindex>=0 && baseindex<fLength) {
  380.         b= (unsigned char) fMasks[baseindex];
  381.         switch (masklevel) {
  382.             case 0: b &= 0x7f; break; // return full mask - top (0x80) bit
  383.             case 1: b &= 1;  break;  
  384.             case 2: b &= 2;  break;
  385.             case 3: b &= 4;  break;
  386.             case 4: b &= 8;  break;
  387.                 // 5-7 reserved for now
  388.             case 5: b &= 16;  break;
  389.             case 6: b &= 32;  break;
  390.             case 7: b &= 64;  break;
  391.             default: b = kMaskEmpty; break;
  392.             }
  393.         }
  394.     return b;
  395. }
  396.  
  397. Boolean DSequence::IsMasked(unsigned char maskbyte, short masklevel)
  398. {
  399.     short b= (unsigned char) maskbyte;
  400.     switch (masklevel) {
  401.         case 0: b &= 0x7f; break; // return full mask - top (0x80) bit
  402.         case 1: b &= 1;  break;  
  403.         case 2: b &= 2;  break;
  404.         case 3: b &= 4;  break;
  405.         case 4: b &= 8;  break;
  406.                 // 5-7 reserved for now
  407.         case 5: b &= 16;  break;
  408.         case 6: b &= 32;  break;
  409.         case 7: b &= 64;  break;
  410.         default: b = kMaskEmpty; break;
  411.         }
  412.     return (b > 0);
  413. }
  414.  
  415.  
  416. inline void SetMaskByte( char& b, short masklevel, short maskval)
  417. {
  418.     switch (masklevel) {
  419.         case 1: if (maskval) b |= 1; else b &= ~1; break;
  420.         case 2: if (maskval) b |= 2; else b &= ~2; break;
  421.         case 3: if (maskval) b |= 4; else b &= ~4; break;
  422.         case 4: if (maskval) b |= 8; else b &= ~8; break;
  423.             // reserve bits 5, 6, 7 & 8 for now 
  424.         case 5: if (maskval) b |= 16; else b &= ~16; break;
  425.         case 6: if (maskval) b |= 32; else b &= ~32; break;
  426.         case 7: if (maskval) b |= 64; else b &= ~64; break;
  427.         // -- 7 & 8 may go unused to store data in printchar form
  428.         default: break;
  429.         }
  430.     //if (b) masklevel--; // debugging b val
  431. }
  432.  
  433. void DSequence::SetMaskAt(long baseindex, short masklevel, short maskval)
  434. {
  435.     if (fMasks && fMasksOk && baseindex>=0 && baseindex<fLength)  
  436.         ::SetMaskByte(fMasks[baseindex],masklevel, maskval);
  437. }
  438.  
  439.  
  440. void DSequence::SetMask(short masklevel, short maskval)
  441. {
  442.     long  baseindex;
  443.     if (fMasks && fMasksOk) 
  444.      for (baseindex=0; baseindex<fLength; baseindex++) 
  445.         ::SetMaskByte(fMasks[baseindex], masklevel, maskval);
  446. }
  447.  
  448.  
  449. inline void FlipMaskByte( char& b, short masklevel)
  450. {
  451.     switch (masklevel) {
  452.         case 1: b ^= 1;  break;  
  453.         case 2: b ^= 2;  break;
  454.         case 3: b ^= 4;  break;
  455.         case 4: b ^= 8;  break;
  456.                 // 5-7 reserved for now?
  457.         case 5: b ^= 16;  break;
  458.         case 6: b ^= 32;  break;
  459.         case 7: b ^= 64;  break;
  460.         default:  break;
  461.         }
  462.     //if (b) masklevel--; // debugging b val
  463. }
  464.  
  465. void DSequence::FlipMaskAt(long baseindex, short masklevel)
  466. {
  467.     if (fMasks && fMasksOk && baseindex>=0 && baseindex<fLength)  
  468.         ::FlipMaskByte(fMasks[baseindex],masklevel);
  469. }
  470.  
  471. void DSequence::FlipMask(short masklevel)
  472. {
  473.     long  baseindex;
  474.     if (fMasks && fMasksOk) 
  475.      for (baseindex=0; baseindex<fLength; baseindex++)  
  476.         FlipMaskByte(fMasks[baseindex],masklevel);
  477. }
  478.  
  479.  
  480. void DSequence::SetInfo( char* theInfo)
  481. {
  482.     if (fInfo) MemFree(fInfo);
  483.     fInfo = StrDup(theInfo);  
  484. }
  485.  
  486. void DSequence::SetName( char* theName)
  487. {
  488.     if (fName) MemFree(fName);
  489.     fName = StrDup(theName);  
  490. }
  491.  
  492. void DSequence::SetKind( short theKind)
  493. {
  494.     fKind= theKind;  
  495. }
  496.  
  497. void DSequence::SetIndex( short theIndex)
  498. {
  499.     fIndex= theIndex;  
  500. }
  501.  
  502. char* DSequence::KindStr() const 
  503. {
  504.     switch (fKind) {
  505.         case kDNA             : return "DNA";
  506.         case kRNA             : return "RNA";
  507.         case kNucleic   : return "Nucleic";
  508.         case kAmino         : return "Protein";
  509.         default                    :
  510.         case kOtherSeq  : return "Unknown";
  511.         }
  512. }
  513.  
  514.  
  515. void DSequence::SetTime( unsigned long theTime)
  516. {
  517.     fModTime= theTime;
  518. }
  519.  
  520. void DSequence::SetOrigin( long theOrigin)
  521. {
  522.     fOrigin= theOrigin;
  523. }
  524.  
  525. void DSequence::Modified()  
  526. {
  527.     fModTime    = gUtil->time();
  528.     fChanged    = TRUE;
  529. }
  530.  
  531. void DSequence::Reformat(short format)
  532. {
  533.     // later...? 
  534. }
  535.  
  536. void DSequence::UpdateFlds()
  537. {
  538.     long         seqlen;
  539.     short     kind;
  540.     unsigned long check;
  541.     fLength= StrLen(fBases);
  542.     GetSeqStats(fBases, fLength, &kind, &check, &seqlen);
  543.     //fValidLength= seqlen; //?? is fLength == StrLen(fBases) or is it <= that, only valid bases??
  544.     fKind            = kind;
  545.     fChecksum    = check;
  546.     //fChanged    = FALSE;
  547. }
  548.  
  549.  
  550.  
  551.  
  552.  
  553. void DSequence::SetSelection( long start, long nbases)
  554. {
  555.     fSelStart= start;
  556.     fSelBases= nbases;
  557. }
  558.  
  559. void DSequence::ClearSelection()
  560. {
  561.     fSelStart= 0;
  562.     fSelBases= 0;
  563. }
  564.  
  565. void    DSequence::GetSelection( long& start, long& nbases)
  566. {
  567.     if (fSelBases==0) {
  568.         start= 0;
  569.         nbases= fLength;
  570.         }    
  571.     else if (fLength>0) {
  572.         start= fSelStart;
  573.         nbases= Min( fSelBases, fLength - fSelStart);
  574.         }    
  575.     else {
  576.         start= 0;
  577.         nbases= 0;
  578.         }
  579. }
  580.  
  581.  
  582. DSequence* DSequence::CopyRange()
  583. {  
  584.             //note: start == 0 == 1st position 
  585.     DSequence* aSeq= (DSequence*) Clone();
  586.     if (fSelStart > 0 || fSelBases > 0) {
  587.         char* h= aSeq->fBases;
  588.         char* m= aSeq->fMasks;
  589.         Boolean domask= (fMasks && fMasksOk);
  590.  
  591.         long len= aSeq->fLength;
  592.         if (fSelStart > 0) {
  593.             len -= fSelStart;
  594.             MemMove( h, h + fSelStart, len);
  595.             if (domask) MemMove( m, m + fSelStart, len);
  596.             }
  597.         if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  598.         h= (char*) MemMore( h, len+1);
  599.         h[len]= '\0'; //((char*)h+len) = '\0';
  600.         aSeq->fBases= h;
  601.         if (domask) {
  602.             m= (char*) MemMore( m, len+1);
  603.             m[len]= '\0';    
  604.             aSeq->fMasks= m;
  605.             }
  606.         }
  607.     aSeq->UpdateFlds();
  608.     aSeq->Modified();
  609.     return aSeq;
  610. }
  611.  
  612.  
  613. DSequence* DSequence::CutRange()
  614. {
  615. /* aSeq.Cutrange means:  
  616.         (a) selected range (fSelStart..fSelStart+fSelBases) is copied into Result; &&
  617.         (b) selected range is removed from this
  618. */
  619.  
  620.     DSequence* aSeq= CopyRange();
  621.     if (fSelStart > 0 || fSelBases > 0) {
  622.         if (fSelStart<0) fSelStart= 0;
  623.         char* h= this->fBases;
  624.         char* m= this->fMasks;
  625.         Boolean domask= (fMasks && fMasksOk);
  626.         long len= this->fLength;
  627.         long top= len - fSelStart - fSelBases;
  628.         if (fSelBases > 0 && fSelBases < len) {
  629.             MemMove( h+fSelStart, h+fSelStart+fSelBases, top);
  630.             if (domask) MemMove( m+fSelStart, m+fSelStart+fSelBases, top);
  631.             }
  632.         len= fSelStart + top;
  633.         h= (char*) MemMore( h, len+1);
  634.         *((char*)h+len)= '\0';       //h[len]= '\0';
  635.         this->fBases= h;
  636.         if (domask) {
  637.             m= (char*) MemMore( m, len+1);
  638.             *((char*)m+len)= '\0';    
  639.             this->fMasks= m;
  640.             }
  641.         }
  642.     UpdateFlds();    
  643.     Modified();
  644.     return aSeq;
  645. }
  646.  
  647. void DSequence::InsertSpacers( long insertPoint, long insertCount, char spacer)
  648. { 
  649.     long i;
  650.     char *cp;
  651.     char* h= this->fBases;
  652.     char* m= this->fMasks;
  653.     Boolean domask= (fMasks && fMasksOk);
  654.     long len= this->fLength;
  655.     if (insertPoint<0) insertPoint= 0;
  656.     else if (insertPoint>len-1) insertPoint= len-1;
  657.     h= (char*) MemMore( h, len+insertCount+1);
  658.     cp = h+insertPoint;
  659.     MemMove( cp+insertCount, cp, len-insertPoint);
  660.     for ( i= 0; i<insertCount; i++)  *cp++= spacer;
  661.     h[len+insertCount]= '\0';
  662.     this->fBases= h;
  663.     if (domask) {
  664.         m= (char*) MemMore( m, len+insertCount+1);
  665.         cp = m+insertPoint;
  666.         char cmask= *cp;
  667.         MemMove( cp+insertCount, cp, len-insertPoint);
  668.         for ( i= 0; i<insertCount; i++)  *cp++= cmask;
  669.         m[len+insertCount]= '\0';
  670.         this->fMasks= m;
  671.         }
  672.     this->UpdateFlds();    
  673.     this->Modified();
  674. }
  675.  
  676.  
  677. DSequence* DSequence::PasteBases( long insertPoint, char* fromBases, char* fromMasks)
  678. { 
  679.     long add;
  680.     DSequence* aSeq= (DSequence*) Clone();
  681.     if (fromBases) add= StrLen( fromBases);  else add= 0;
  682.     if (add>0) {
  683.         char* h= aSeq->fBases;
  684.         char* m= aSeq->fMasks;
  685.         Boolean domask= (fMasks && fMasksOk);
  686.         long len = aSeq->fLength;
  687.         if (insertPoint<=0) {
  688.             fromBases= StrDup( fromBases);
  689.             StrExtendCat( &fromBases, h);
  690.             MemFree( h);
  691.             h= fromBases;
  692.             if (domask && fromMasks) {
  693.                 fromMasks= StrDup( fromMasks);
  694.                 StrExtendCat( &fromMasks, m);
  695.                 MemFree( m);
  696.                 m= fromMasks;
  697.                 }
  698.             }        
  699.         else if (insertPoint>=len) { 
  700.             StrExtendCat( &h, fromBases);
  701.             if (domask) StrExtendCat( &m, fromMasks);
  702.             }
  703.         else { 
  704.             h= (char*) MemMore( h, len+add);
  705.             MemMove( h+insertPoint+add, h+insertPoint, len-insertPoint);
  706.             MemMove( h+insertPoint+add, fromBases, add);
  707.             *((char*)h+len+add)= '\0';
  708.             if (domask) {
  709.                 m= (char*) MemMore( m, len+add);
  710.                 MemMove( m+insertPoint+add, m+insertPoint, len-insertPoint);
  711.                 MemMove( m+insertPoint+add, fromMasks, add);
  712.                 *((char*)m+len+add)= '\0';
  713.                 }
  714.             }
  715.         aSeq->fBases= h;
  716.         aSeq->fMasks= m;
  717.         }
  718.     aSeq->UpdateFlds();    
  719.     aSeq->Modified();
  720.     return aSeq;
  721. }
  722.  
  723.  
  724. DSequence* DSequence::Reverse()
  725. {
  726.     long i;
  727.     DSequence* aSeq= (DSequence*) Clone();    
  728.     if (fSelStart < 0) fSelStart= 0;
  729.     long len= aSeq->fLength - fSelStart;
  730.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  731.     char* pc= this->fBases+fSelStart + len;
  732.     char* pr= aSeq->fBases+fSelStart;
  733.     for (i= 0; i<len; i++) *pr++= *--pc;
  734.     Boolean domask= (fMasks && fMasksOk);
  735.     if (domask) {
  736.         pc= this->fMasks+fSelStart + len;
  737.         pr= aSeq->fMasks+fSelStart;
  738.         for (i= 0; i<len; i++) *pr++= *--pc;
  739.         }
  740.     aSeq->UpdateFlds();
  741.     aSeq->Modified();
  742.     return aSeq;
  743. }
  744.  
  745.  
  746. DSequence* DSequence::Complement()
  747. // A->T, T->A, G->C, C->G  
  748. {
  749.     Boolean        isrna;
  750.     DSequence* aSeq= (DSequence*) Clone();    
  751.     
  752.     if (fKind == kRNA) isrna= TRUE;
  753.     else if (fKind == kDNA || fKind == kNucleic) isrna= FALSE;
  754.     else return aSeq;    
  755.         
  756.     if (fSelStart < 0) fSelStart= 0;
  757.     long len= aSeq->fLength - fSelStart;
  758.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  759.     char* pc= this->fBases+fSelStart;
  760.     char* pr= aSeq->fBases+fSelStart;
  761.     NucleicComplement( isrna, pc, pr, len);
  762.         
  763.     aSeq->UpdateFlds();
  764.     aSeq->Modified();
  765.     return aSeq;
  766. }
  767.  
  768.  
  769.  
  770. DSequence* DSequence::ChangeCase(Boolean toUpper)
  771. // g->G or G->g
  772. {
  773.     DSequence* aSeq= (DSequence*) Clone();    
  774.     
  775.     if (fSelStart < 0) fSelStart= 0;
  776.     char* h= aSeq->fBases;
  777.     long len= aSeq->fLength - fSelStart;
  778.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  779.     char* pc= this->fBases+fSelStart;
  780.     char* pr= aSeq->fBases+fSelStart;
  781.     
  782.     long i;
  783.     if (toUpper) for (i= 0; i<len; i++) {
  784.         *pr++= toupper(*pc); pc++;
  785.         }
  786.     else for (i= 0; i<len; i++) {
  787.         *pr++= tolower(*pc); pc++;
  788.         }    
  789.         
  790.     aSeq->UpdateFlds();
  791.     aSeq->Modified();
  792.     return aSeq;
  793. }
  794.  
  795. DSequence* DSequence::Dna2Rna(Boolean toRna)
  796. // T->U or U->T 
  797. {
  798.     Boolean        isrna;
  799.     DSequence* aSeq= (DSequence*) Clone();    
  800.     
  801.     if (fKind==kRNA) isrna= TRUE;
  802.     else if (fKind == kDNA || fKind == kNucleic) isrna= FALSE;
  803.     else return aSeq;    
  804.      
  805.     if (fSelStart < 0) fSelStart= 0;
  806.     long len= aSeq->fLength - fSelStart;
  807.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  808.     char* pc= this->fBases+fSelStart;
  809.     char* pr= aSeq->fBases+fSelStart;
  810.     
  811.     long i;
  812.     char b, c;
  813.     if (toRna) for (i= 0; i<len; i++) {
  814.         b= *pc++;
  815.         if (b=='T') c= 'U';
  816.         else if (b=='t') c= 'u';
  817.         else c= b;
  818.         *pr++= c;
  819.         }
  820.     else for (i= 0; i<len; i++) {
  821.         b= *pc++;
  822.         if (b=='U') c= 'T';
  823.         else if (b=='u') c= 't';
  824.         else c= b;
  825.         *pr++= c;
  826.         }    
  827.         
  828.     aSeq->UpdateFlds();
  829.     aSeq->Modified();
  830.     return aSeq;
  831. }
  832.  
  833.  
  834. DSequence* DSequence::LockIndels( Boolean doLock)
  835. // indelSoft -> indelHard or vice versa 
  836. {    
  837.     char b, c, fromc, toc;
  838.     DSequence* aSeq= (DSequence*) Clone();    
  839.     
  840.     if (fSelStart < 0) fSelStart= 0;
  841.     long  len= aSeq->fLength - fSelStart;
  842.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  843.     char* pc= this->fBases+fSelStart;
  844.     char* pr= aSeq->fBases+fSelStart;
  845.     
  846.     if (doLock) { fromc= indelSoft; toc= indelHard; }
  847.     else { fromc= indelHard; toc= indelSoft; }
  848.     
  849.     for (long i= 0; i<len; i++) {
  850.         b= *pc++;
  851.         if (b == fromc) c= toc;
  852.         else c= b;
  853.         *pr++= c;
  854.         }    
  855.         
  856.     aSeq->UpdateFlds();
  857.     aSeq->Modified();
  858.     return aSeq;
  859. }
  860.  
  861.  
  862. DSequence* DSequence::Slide( long slideDist)
  863. /* insert "-" below (+dist) or above (-dist) of seq range,
  864.     squeeze out "-" above (below) of range
  865. */
  866. {     
  867.     long         slideSave, newlen, fulllen, len, i, j= 0;
  868.     char         b, cutchar, cutmask;
  869.     char        *pc, *pr;
  870.     
  871.     DSequence* aSeq= (DSequence*) Clone();    
  872.     if (fSelStart < 0) fSelStart= 0;
  873.     fulllen= aSeq->fLength;
  874.     len= fulllen - fSelStart;
  875.     Boolean domask= (fMasks && fMasksOk);
  876.     Boolean dosqueeze= true; //(indelSoft & 0x80 != 0); // quick fix flag
  877.     indelSoft &= 0x7f;
  878.     
  879.     if (slideDist < 0) {
  880.         // tough... do if below works...
  881.         // really need ability to extend seq in 5' (lower) direction 
  882.         //! build pr in reverse, then flip it....
  883.  
  884.         if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  885.         slideDist= -slideDist;
  886.         slideSave= slideDist;
  887.         aSeq->fBases= (char*) MemMore(aSeq->fBases, 2+fulllen+slideDist);
  888.         pr= aSeq->fBases;
  889.         pc= this->fBases + fulllen;
  890.         for (i= fulllen-1; i >= fSelStart+len; i--) *pr++= *--pc;
  891.              
  892.         if (pr == aSeq->fBases || *(pr-1) == indelEdge) cutchar= indelEdge;
  893.         else cutchar= indelSoft; //??
  894.         for (i=0; i<slideDist; i++) *pr++= cutchar;
  895.  
  896.         pc= this->fBases+fSelStart+len;
  897.         for (i= fSelStart+len-1; i >= fSelStart; i--) *pr++= *--pc;
  898.             
  899.         pc= this->fBases+fSelStart;
  900.         for (i= fSelStart-1; i>=0; i--) {
  901.             b= *--pc;
  902.             if (dosqueeze && slideDist>0 && (b==indelSoft || b==indelEdge)) slideDist--;
  903.             else {
  904.                 *pr++= b;
  905.                 }
  906.             }
  907.  
  908.         newlen= pr - aSeq->fBases;
  909.         char* htmp= (char*) MemNew(newlen+1);
  910.         pc= htmp;
  911.         for (i=0; i<newlen; i++) *pc++= *--pr; //pr[j-i+1];
  912.         pc= '\0';
  913.         MemFree(aSeq->fBases);
  914.         aSeq->fBases = htmp;  
  915.         aSeq->fLength= newlen;
  916.         aSeq->fOrigin= fulllen - newlen; //== this->fLength - aSeq->fLength
  917.  
  918.         if (domask) {
  919.             slideDist= slideSave;
  920.             aSeq->fMasks= (char*) MemMore(aSeq->fMasks, 2+fulllen+slideDist);
  921.             pr= aSeq->fMasks;
  922.             pc= this->fMasks + fulllen;
  923.             for (i= fulllen-1; i >= fSelStart+len; i--) *pr++= *--pc;
  924.                  
  925.             cutmask= 0x80; // required empty mask  
  926.             for (i=0; i<slideDist; i++) *pr++= cutmask;
  927.     
  928.             pc= this->fMasks+fSelStart+len;
  929. #ifdef WIN_MSWIN
  930.             // mswin is crashing here for no obvious reason
  931.             //pc--; len--;
  932.             for (i= fSelStart + len - 1; i >= fSelStart; i--) {
  933.                  register short amask;
  934.                  --pc;
  935.                  amask= *pc;    //mswin bomb is here !
  936.                  *pr= amask;
  937.                  pr++;
  938.          }
  939. #else
  940.             for (i= fSelStart+len-1; i >= fSelStart; i--) *pr++= *--pc;
  941. #endif
  942.             pc= this->fMasks+fSelStart;
  943.             for (i= fSelStart-1; i>=0; i--) {
  944.                 b= *--pc;
  945.                 if (slideDist>0 && (b==cutmask)) slideDist--;
  946.                 else {
  947.                     *pr++= b;
  948.                     }
  949.                 }
  950.                 
  951.             newlen= pr - aSeq->fMasks;
  952.             htmp= (char*) MemNew(newlen+1);
  953.             pc= htmp;
  954.             for (i=0; i<newlen; i++) *pc++= *--pr; //pr[j-i+1];
  955.             pc= '\0';
  956.             MemFree(aSeq->fMasks);
  957.             aSeq->fMasks = htmp;  
  958.             }
  959.             
  960.         }    
  961.         
  962.     else if (slideDist > 0) {
  963.         fulllen= len;
  964.         if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  965.         aSeq->fBases= (char*) MemMore(aSeq->fBases, 2+fSelStart+fulllen+slideDist); //make it large enough
  966.  
  967.         pr= aSeq->fBases;
  968.       if (fSelStart==0 || pr[fSelStart-1]==indelEdge) cutchar= indelEdge; 
  969.         else cutchar= indelSoft; 
  970.  
  971.         pc= this->fBases+fSelStart;
  972.         pr= aSeq->fBases+fSelStart;
  973.  
  974.         for (j= 0; j<slideDist; j++) *pr++= cutchar;
  975.         for (i= 0; i<len; i++) *pr++= *pc++; 
  976.     
  977.         for (i= len; i<fulllen; i++) {
  978.             b= *pc++;
  979.             if (dosqueeze && slideDist>0 && (b==indelSoft || b==indelEdge)) 
  980.                 slideDist--;
  981.             else {
  982.                 *pr++= b;
  983.                 }
  984.             }
  985.         
  986.         *pr= 0;
  987.         newlen= pr - aSeq->fBases; 
  988.         aSeq->fBases= (char*)MemMore(aSeq->fBases, newlen+1);  
  989.         
  990.         if (domask) {
  991.             aSeq->fMasks= (char*) MemMore(aSeq->fMasks, 2+fSelStart+fulllen+slideDist);
  992.           cutmask= 0x80; 
  993.             pc= this->fMasks+fSelStart;
  994.             pr= aSeq->fMasks+fSelStart;
  995.     
  996.             for (j= 0; j<slideDist; j++) *pr++= cutmask;
  997.             for (i= 0; i<len; i++) *pr++= *pc++; 
  998.  
  999.             for (i= len; i<fulllen; i++) {
  1000.                 b= *pc++;
  1001.                 if (slideDist>0 && (b==cutmask)) // this is BAD? must be done w/ base slide !?
  1002.                     slideDist--;
  1003.                 else {
  1004.                     *pr++= b;
  1005.                     }
  1006.                 }
  1007.             
  1008.             *pr= 0;
  1009.             newlen= pr - aSeq->fMasks; 
  1010.             aSeq->fMasks= (char*)MemMore(aSeq->fMasks, newlen+1);  
  1011.             }
  1012.         }
  1013.             
  1014.     aSeq->UpdateFlds();
  1015.     aSeq->Modified();
  1016.     return aSeq;
  1017. }
  1018.  
  1019.  
  1020. DSequence* DSequence::Compress()
  1021. // pull out "-", ".", ? anything but nuc/aa codes ? 
  1022. {
  1023.     long  len, fulllen, i;
  1024.     char     b, c;
  1025.     char * mc, * mr;
  1026.     Boolean domask= (fMasks && fMasksOk);
  1027.  
  1028.     DSequence* aSeq= (DSequence*) Clone();    
  1029.     if (fSelStart < 0) fSelStart= 0;
  1030.     len= aSeq->fLength - fSelStart;
  1031.     fulllen= len;
  1032.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  1033.     char* pc= this->fBases+fSelStart;
  1034.     char* pr= aSeq->fBases+fSelStart;
  1035.     if (domask) {
  1036.         mc= this->fMasks+fSelStart;
  1037.         mr= aSeq->fMasks+fSelStart;
  1038.         }
  1039.         
  1040.     if (domask) for (i= 0; i<len; i++) {
  1041.         b= *pc++;
  1042.         c= *mc++;
  1043.         if (b!=indelSoft && b!=indelEdge) { 
  1044.             *pr++= b;
  1045.             *mr++= c;
  1046.             }
  1047.         }
  1048.     else for (i= 0; i<len; i++) {
  1049.         b= *pc++;
  1050.         if (b!=indelSoft && b!=indelEdge) 
  1051.             *pr++= b;
  1052.         }
  1053.     for (i=len; i<fulllen; i++) *pr++= *pc++;
  1054.     if (domask) for (i=len; i<fulllen; i++) *mr++= *mc++;
  1055.  
  1056.     long newlen= pr - aSeq->fBases;
  1057.     aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen+1);
  1058.     aSeq->fBases[newlen]= 0;
  1059.     if (domask) {
  1060.         aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen+1);
  1061.         aSeq->fMasks[newlen]= 0;
  1062.         }
  1063.         
  1064.     aSeq->UpdateFlds();
  1065.     aSeq->Modified();
  1066.     return aSeq;
  1067. }
  1068.  
  1069.  
  1070.  
  1071.  
  1072. DSequence* DSequence::CompressFromMask(short masklevel)
  1073. {
  1074.     long  len, fulllen, i;
  1075.     char     b, c;
  1076.     char * mc, * mr;
  1077.     Boolean domask= (fMasks && fMasksOk && masklevel>0);
  1078.  
  1079.     if (!domask) return NULL; // or aSeq ???
  1080.     DSequence* aSeq= (DSequence*) Clone();    
  1081.     //if (!domask) return aSeq;  
  1082.  
  1083.     if (fSelStart < 0) fSelStart= 0;
  1084.     len= aSeq->fLength - fSelStart;
  1085.     fulllen= len;
  1086.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  1087.     char* pc= this->fBases+fSelStart;
  1088.     char* pr= aSeq->fBases+fSelStart;
  1089.     mc= this->fMasks+fSelStart;
  1090.     mr= aSeq->fMasks+fSelStart;
  1091.          
  1092.     for (i= 0; i<len; i++) {
  1093.         b= *pc++;
  1094.         c= *mc++;
  1095.         if (IsMasked(c, masklevel)) { 
  1096.             *pr++= b;
  1097.             *mr++= c;
  1098.             }
  1099.         }
  1100.     for (i=len; i<fulllen; i++) *pr++= *pc++;
  1101.     for (i=len; i<fulllen; i++) *mr++= *mc++;
  1102.  
  1103.     long newlen= pr - aSeq->fBases;
  1104.     aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen+1);
  1105.     aSeq->fBases[newlen]= 0;
  1106.     aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen+1);
  1107.     aSeq->fMasks[newlen]= 0;
  1108.         
  1109.     aSeq->UpdateFlds();
  1110.     aSeq->Modified();
  1111.     return aSeq;
  1112. }
  1113.  
  1114.  
  1115.  
  1116. DSequence* DSequence::Translate(Boolean keepUnselected)
  1117. {    
  1118.     char  c;
  1119.     char    acodon[3];
  1120.     Boolean isrna, nocodon;
  1121.     char    *pr, *pc, *mr, *mc;
  1122.     long    len, fulllen, newlen, i, j, k;
  1123.     
  1124.     DSequence* aSeq= (DSequence*) Clone();        
  1125.     if (DCodons::NotAvailable()) return aSeq; // safer than returning NULL !?
  1126.     // !? what do we do w/ fMasks !?
  1127.     Boolean domask= (fMasks && fMasksOk);
  1128.     
  1129.     if (fSelStart < 0) fSelStart= 0;
  1130.     fulllen= this->fLength - fSelStart;
  1131.     len= fulllen;
  1132.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  1133.     
  1134.     switch (fKind) {
  1135.     
  1136.         case kRNA:
  1137.         case kDNA:
  1138.         case kNucleic: 
  1139.             {
  1140.             if (keepUnselected) {
  1141.                 pr= aSeq->fBases + fSelStart; // start at selection, leaving untranslated 5'
  1142.                 if (domask) mr= aSeq->fMasks + fSelStart;
  1143.                 }
  1144.             else {
  1145.                 pr= aSeq->fBases;
  1146.                  if (domask) mr= mr= aSeq->fMasks;
  1147.                  }
  1148.                  
  1149.             pc= this->fBases + fSelStart;
  1150.             if (domask) mc= this->fMasks + fSelStart;
  1151.             isrna= (fKind == kRNA);
  1152.             i= fSelStart;
  1153.             while (i < len) {
  1154.                 for (j=0; j<3; j++) {
  1155.                     c= toupper( *pc); pc++; // toupper is sometimes a MACRO ! no (*pc++) !
  1156.                     if (isrna && c=='U') c= 'T';
  1157.                     acodon[j]= c;
  1158.                     }
  1159.                 i += 3;
  1160.                 for (k= 0, nocodon= true; k<64; k++) 
  1161.                     if ( MemCmp(DCodons::codon(k),acodon,3) == 0) { 
  1162.                         *pr++= DCodons::amino(k); //gCodonTable[k].amino;
  1163.                         nocodon= false;
  1164.                         break;
  1165.                         }
  1166.                 if (nocodon) *pr++= '?';
  1167.                 if (domask) {
  1168.                     *mr++= *mc;
  1169.                     mc += 3;
  1170.                     }
  1171.                 }
  1172.                  // copy over untranslated portion
  1173.             if (keepUnselected) {
  1174.                 long savei= i;
  1175.                 for ( ; i < fulllen; i++)  *pr++= *pc++; 
  1176.                 if (domask) for ( i= savei; i < fulllen; i++)  *mr++= *mc++; 
  1177.                 }
  1178.             }
  1179.             break;
  1180.             
  1181.         case kAmino: 
  1182.             {
  1183.             char *bestcodon;
  1184.             long  selend;
  1185.             
  1186.             if (len < fulllen) {
  1187.                 newlen= fulllen + len*3;
  1188.                 selend= fSelStart + len;
  1189.                 }
  1190.             else {
  1191.                 newlen= fulllen * 3;
  1192.                 selend= len;
  1193.                 }
  1194.             pc= this->fBases + fSelStart;
  1195.             aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen);
  1196.             if (keepUnselected)
  1197.                 pr= aSeq->fBases + fSelStart;
  1198.             else
  1199.                 pr= aSeq->fBases;
  1200.                 
  1201.             if (domask) {
  1202.                 mc= this->fMasks + fSelStart;
  1203.                 aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen);
  1204.                 if (keepUnselected)
  1205.                     mr= aSeq->fMasks + fSelStart;
  1206.                 else
  1207.                     mr= aSeq->fMasks;
  1208.                 }
  1209.                 
  1210.             for (i= fSelStart; i<selend; i++) {
  1211.                 c= toupper( *pc); pc++;// toupper is sometimes a MACRO ! no (*pc++) !
  1212.                 bestcodon= (char*) DCodons::FindBestCodon(c);
  1213.                 MemCpy( pr, bestcodon, 3);  
  1214.                 pr += 3;
  1215.                 if (domask) {
  1216.                     for (j=0; j<3; j++) *mr++= *mc;
  1217.                     mc++;
  1218.                     }
  1219.                 }
  1220.                  // copy over untranslated portion
  1221.             if (keepUnselected) {
  1222.                 long savei= i;
  1223.                 for ( ; i < fulllen; i++)  *pr++= *pc++; 
  1224.                 if (domask) for ( i=savei; i < fulllen; i++)  *mr++= *mc++; 
  1225.                 }
  1226.             }
  1227.             break;
  1228.         
  1229.         default:
  1230.             return NULL;
  1231.          
  1232.       }
  1233.         
  1234.     *pr= '\0';
  1235.     newlen= pr - aSeq->fBases;
  1236.     aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen);
  1237.     if (domask) {
  1238.         *mr= '\0';
  1239.         newlen= mr - aSeq->fMasks;
  1240.         aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen);
  1241.         }
  1242.     aSeq->UpdateFlds();
  1243.     aSeq->Modified();
  1244.     return aSeq;    
  1245. }
  1246.  
  1247.  
  1248. DSequence* DSequence::OnlyORF(char nonorf)
  1249. // convert sequence to ORF part only 
  1250. {
  1251.     DSequence* aSeq= (DSequence*) Clone();    
  1252.     
  1253.     if (fSelStart < 0) fSelStart= 0;
  1254.     long len= aSeq->fLength - fSelStart;
  1255.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  1256.  
  1257.     long start = 0, fini= 0;
  1258.     long anchor= fSelStart;
  1259.     do {
  1260.         aSeq->SearchORF( start, fini);
  1261.         if (start != kSearchNotFound) {
  1262.             char* pr= aSeq->fBases + anchor;
  1263.             if (start>len) start= len;
  1264.             if (fini<start) fini= start+1;
  1265.             for ( ; anchor<start; anchor++) *pr++ = nonorf; 
  1266.             anchor= fini+1;
  1267.             fSelStart= anchor;
  1268.             }
  1269.     } while (start != kSearchNotFound);
  1270.         
  1271.     aSeq->UpdateFlds();
  1272.     aSeq->Modified();
  1273.     return aSeq;
  1274. }
  1275.  
  1276.  
  1277. void DSequence::SetORFmask(short masklevel, short frame)
  1278. {
  1279.     if (!fMasks || !fMasksOk) {
  1280.         FixMasks();
  1281.         if (!fMasks || !fMasksOk) return;
  1282.         }
  1283.         
  1284.     // revise this to do multiple frames in one go,
  1285.     //  offseting each anchor +1 beyond START of last orf, not after STOP
  1286.     // since separate frame calls return the same orfs !
  1287.     
  1288.     frame %= 3; // drop reverse frames
  1289.     long savesel= fSelStart;
  1290.     long savelen= fSelBases;
  1291.     fSelStart= frame;
  1292.     fSelBases= fLength - frame;
  1293.  
  1294.     long start = 0, fini= 0;
  1295.     long anchor= fSelStart;
  1296.     do {
  1297.         this->SearchORF( start, fini);
  1298.         if (start != kSearchNotFound) {
  1299.             long i;
  1300.             for (i= anchor; i<start; i++) SetMaskAt( i, masklevel, 0);
  1301.             for (i= start; i<=fini ; i++) SetMaskAt( i, masklevel, 1);
  1302.             if (fini == kSearchNotFound) {
  1303.                 start= kSearchNotFound;
  1304.                 }
  1305.             else {
  1306.               anchor= fini+1;
  1307.               fSelStart= anchor; // start+1 // !? want overlapped orfs !?
  1308.               }
  1309.             }
  1310.     } while (start != kSearchNotFound);
  1311.     
  1312.     fSelStart= savesel;
  1313.     fSelBases= savelen;
  1314. }
  1315.  
  1316. void DSequence::SetShowORF( short turnon) 
  1317. { 
  1318.     fShowORF= turnon > 0; 
  1319.     if (turnon > 1) {
  1320.         if (turnon & 2) SetORFmask(5, 0);
  1321.         //if (turnon & 4) SetORFmask(6, 1);
  1322.         //if (turnon & 8) SetORFmask(7, 2);
  1323.         }
  1324. }
  1325.  
  1326. void DSequence::SetShowRE( short turnon) 
  1327. { fShowRE= turnon > 0; 
  1328. }
  1329.  
  1330. void DSequence::SetShowTrace( short turnon) 
  1331. { fShowTrace= turnon > 0; 
  1332. }
  1333.  
  1334.  
  1335. //Boolean domask= (fMasks && fMasksOk);
  1336.  
  1337.  
  1338. void DSequence::DoWrite(DFile* aFile, short format)
  1339. {  
  1340.     if (fLength>0) {
  1341.         DSeqFile::WriteSeqWrapper( aFile, fBases, fLength, format, fName);
  1342.         if (fMasks && fMasksOk && DSeqFile::gWriteMasks) 
  1343.             DSeqFile::WriteMaskWrapper( aFile, fMasks, fLength, format, fName);
  1344.         }
  1345. }
  1346.  
  1347.  
  1348. void DSequence::DoWriteSelection(DFile* aFile, short format)
  1349. {
  1350.     if (fSelBases==0) 
  1351.         DoWrite( aFile, format);
  1352.     else if (fLength>0) {
  1353.         long aStart= Min( fSelStart, fLength);
  1354.         long aLength= Min( fSelBases, fLength - fSelStart);
  1355.         DSeqFile::WriteSeqWrapper(  aFile, fBases+aStart, aLength, format, fName);
  1356.         if (fMasks && fMasksOk && DSeqFile::gWriteMasks)  
  1357.             DSeqFile::WriteMaskWrapper( aFile, fMasks+aStart, aLength, format, fName);
  1358.         }
  1359. }
  1360.  
  1361.  
  1362.  
  1363.  
  1364. Boolean DSequence::GoodChar(char& aChar) 
  1365. {
  1366.     char* aacodes = gAminos;
  1367.     char* iupaccodes = gIubbase;
  1368.     char* nacodes = gPrimenuc;
  1369.     char    c;
  1370.     char* codes = NULL;
  1371.  
  1372.   if (aChar=='\n' || aChar=='\r')  
  1373.         return FALSE;
  1374.     else if (aChar<' ' || aChar>'~')  
  1375.         return TRUE; //allow return,tab,delete, ? option chars
  1376.  
  1377.   switch (fKind) {
  1378.       case kDNA            :
  1379.         case kRNA            : codes= nacodes; break;
  1380.         case kNucleic    : codes= iupaccodes; break;
  1381.         case kAmino        : codes= aacodes; break;
  1382.         default:
  1383.       case kOtherSeq: return TRUE;
  1384.         }
  1385.         
  1386.   c= toupper(aChar);
  1387.     while (*codes) {
  1388.         if (c == *codes) {
  1389.           switch (fKind) {
  1390.                 case kDNA: if (c == 'U') {
  1391.                     if (aChar=='u') aChar= 't'; else aChar= 'T';
  1392.                     }
  1393.                 break;
  1394.                 case kRNA: if (c == 'T') {
  1395.                     if (aChar=='t') aChar= 'u'; else aChar= 'U';
  1396.                     }
  1397.                 break;
  1398.                 }
  1399.             return TRUE;
  1400.             }
  1401.         codes++;
  1402.         }
  1403.     return FALSE;
  1404. }
  1405.  
  1406.  
  1407. Boolean DSequence::IsConsensus()
  1408. {
  1409.     return (fName && StringCmp(fName,kConsensus)==0);
  1410. }
  1411.  
  1412.  
  1413.  
  1414.  
  1415.  
  1416. long DSequence::SearchAgain() 
  1417. {
  1418.     return Search( "", FALSE);
  1419. }
  1420.  
  1421.  
  1422.  
  1423. void DSequence::SearchORF( long& start, long& stop) 
  1424. {
  1425.     char  startc[2];
  1426.     long  i, j, k;
  1427.     short    ns, aaMinORFsize;
  1428.     CodonStat  stops[10];
  1429.     
  1430.     start= kSearchNotFound;
  1431.     stop = kSearchNotFound;
  1432.     if (DCodons::NotAvailable()) return;
  1433.     
  1434.     for (ns= 0, i= 0; i<64; i++) {
  1435.         if (ns<10 && DCodons::amino(i) == '*') 
  1436.             stops[ns++] = DCodons::fCodons[i];
  1437.         }
  1438.     if (ns == 0) return;
  1439.         
  1440.     switch (fKind) {  
  1441.         case kRNA:
  1442.         case kDNA:
  1443.         case kNucleic: 
  1444.              start= Search( DCodons::startcodon(), FALSE);
  1445.              if (start == kSearchNotFound) return;
  1446.              for (j= start+3; j<fLength; j += 3) {
  1447.                if ( j+3 - start >= DSeqList::gMinORFsize) 
  1448.                    for (k=0; k<ns; k++) {
  1449. #if 1
  1450.                        short jc, jt;
  1451.                        for (jc= 0, jt=0; jc<3; ) {
  1452.                            long tbit= NucleicBits( fBases[j+jt]);
  1453.                       long sbit= NucleicBits( stops[k].codon[jc]); 
  1454.                       if (tbit == kMaskIndel) jt++;
  1455.                         else if (!(kMaskNucs & tbit & sbit)) goto nextcodon;
  1456.                         else { jt++; jc++; }
  1457.                         }
  1458.                     stop= j+3; 
  1459.                     return;
  1460. #else                    
  1461.                      if (StrNICmp( fBases+j, stops[k].codon, 3) == 0) {
  1462.                          stop= j+3;
  1463.                          return;
  1464.                          }
  1465. #endif                
  1466.          nextcodon:
  1467.                      continue;
  1468.                      }
  1469.                  }
  1470.             break;
  1471.             
  1472.         case kAmino:
  1473.             startc[0]= DCodons::startamino();
  1474.             startc[1]= 0;
  1475.              aaMinORFsize= (DSeqList::gMinORFsize+2) / 3;
  1476.              start= Search( startc, FALSE);
  1477.              if (start == kSearchNotFound) return;
  1478.              for (j= start+1; j<fLength; j++) {
  1479.                if ( j - start >= aaMinORFsize) 
  1480.                   for (k=0; k<ns; k++) {
  1481.                      if (toupper(fBases[j]) == toupper(stops[k].amino)) {
  1482.                          stop= j+1; // +1 to include this codon
  1483.                          return;
  1484.                          }
  1485.                      }
  1486.                  }
  1487.             break;
  1488.             
  1489.         default:
  1490.             break;
  1491.         }
  1492. }
  1493.  
  1494.  
  1495.  
  1496.  
  1497. void DSequence::NucleicSearch(char* source, char* target, SearchRec& aSearchRec)
  1498. {
  1499.     long k, j, tBits, lens, ti;
  1500.     Boolean done;
  1501.  
  1502.                     //find target[ti] where tBits == kMaskA,C,G,T for fastest search
  1503.                     //This part is wasted time for SearchAgain... store tBits/ti in saverec...
  1504.     lens= StrLen(target);
  1505.     for (j= 0, ti= 0, tBits= 0; j<lens && !tBits; j++) {
  1506.         ti= j;
  1507.         tBits= NucleicBits( target[ti]);
  1508.         tBits &= kMaskNucs;  //drop RNA bit
  1509.         }
  1510.  
  1511.     do {
  1512. reloop:
  1513.             aSearchRec.indx += aSearchRec.step; 
  1514.             aSearchRec.lim  -= aSearchRec.step;    
  1515.             
  1516.             if (aSearchRec.step>0) {
  1517.                 for (j= 0; j<aSearchRec.lim; j++)
  1518.                  // if (IsNucleicMatch( tBits, NucleicBits( source[j+aSearchRec.indx]))) 
  1519.                  if ( tBits & NucleicBits( source[j+aSearchRec.indx]) ) 
  1520.                  {
  1521.                     k= j; 
  1522.                     goto nextindx;
  1523.                  }
  1524.                 k= aSearchRec.lim;
  1525.                 }            
  1526.             else {  //step<0 IS BAD .. see UTextDoc version...
  1527.                 for (j= 0; j > -aSearchRec.lim; j--) 
  1528.                   //if (IsNucleicMatch( tBits, NucleicBits( source[j+aSearchRec.indx])))
  1529.                   if ( tBits & NucleicBits( source[j+aSearchRec.indx]) ) 
  1530.                 {
  1531.                     k= j; 
  1532.                     goto nextindx;
  1533.                     }
  1534.                 k= -aSearchRec.lim;
  1535.                 }
  1536.                 
  1537. nextindx:  
  1538.             aSearchRec.indx += k; 
  1539.             aSearchRec.lim  -= k;
  1540.                 //? not enough source left to match full target 
  1541.             if (aSearchRec.step>0) {
  1542.                 if (aSearchRec.lim < lens-ti) aSearchRec.lim= 0;
  1543.                 }            
  1544.             else {
  1545.                 if (aSearchRec.lim < ti) aSearchRec.lim= 0; //? or lim < ti-1 
  1546.                 }
  1547.             done= (aSearchRec.lim==0);
  1548.             
  1549.             if (!done) {
  1550.                 long srcstart= aSearchRec.indx - ti;
  1551.                 long js, jt;
  1552.                 for (js= 0, jt= 0; jt<lens; ) {
  1553. #if 1
  1554.                   long tbit= NucleicBits( target[jt]);
  1555.                   long sbit= NucleicBits( source[js+srcstart]); 
  1556.                   if (tbit == kMaskIndel) jt++;
  1557.                   else if (sbit == kMaskIndel) js++;
  1558.                     else if (!(kMaskNucs & tbit & sbit)) goto reloop;
  1559.                     else { jt++; js++; }
  1560. #else 
  1561.                   if (!(kMaskNucs & NucleicBits( target[jt]) & NucleicBits( source[js+srcstart]))) 
  1562.                       goto reloop;
  1563.                     else { jt++; js++; }
  1564. #endif
  1565.                     }
  1566.                 done= true;
  1567.                 }
  1568.         } while (!done);
  1569. }  
  1570.     
  1571.  
  1572. void DSequence::ProteinSearch(char* source, char* target, SearchRec& aSearchRec)
  1573. {
  1574.     //just toupper match of chars 
  1575.     char      tChar; 
  1576.     Boolean  done;
  1577.     long        lens, k, j;
  1578.  
  1579.     lens= StrLen( target);
  1580.     tChar= toupper( target[0]);
  1581.     
  1582.     do {
  1583. reloop:
  1584.             aSearchRec.indx += aSearchRec.step; 
  1585.             aSearchRec.lim  -= aSearchRec.step;    
  1586.             
  1587.             if (aSearchRec.step>0) {
  1588.                 for (j= 0; j<aSearchRec.lim; j++) 
  1589.                  if (tChar == toupper( source[j+aSearchRec.indx])) {
  1590.                     k= j;
  1591.                     goto nextindx;
  1592.                     }
  1593.                 k= aSearchRec.lim;
  1594.                 }            
  1595.             else {
  1596.                 for (j= 0; j > -aSearchRec.lim; j--)
  1597.                  if (tChar == toupper(source[j+aSearchRec.indx])) {
  1598.                     k= j; 
  1599.                     goto nextindx;
  1600.                     }
  1601.                 k= -aSearchRec.lim;
  1602.                 }
  1603.                 
  1604. nextindx:  
  1605.             aSearchRec.indx += k; 
  1606.             aSearchRec.lim  -= k;
  1607.                 //? not enough source left to match full target 
  1608.             if (aSearchRec.step>0) {
  1609.                 if (aSearchRec.lim < lens-1) aSearchRec.lim= 0;
  1610.                 }            
  1611.             else {
  1612.                 if (aSearchRec.lim < 1) aSearchRec.lim= 0;  
  1613.                 }
  1614.             done= (aSearchRec.lim==0);
  1615.             
  1616.             if (!done) {
  1617.                 for (j= 1; j<lens; j++) {
  1618.                     if (toupper( target[j]) != toupper( source[j+aSearchRec.indx])) 
  1619.                         goto reloop;
  1620.                     }
  1621.                 done= true;
  1622.                 }
  1623.         } while (!done);
  1624. }  
  1625.  
  1626.  
  1627.     
  1628.  
  1629. long DSequence::Search( char* target, Boolean backwards)
  1630. /* search for target in sequence;
  1631.     search == abs index into fBases[0..fLength] or kSearchNotFound
  1632.   if (target='') then find next  
  1633.   if (backwards) then backwards search
  1634. */
  1635. {
  1636.     long        limit, aStart;
  1637.     SearchRec  aSearchRec;
  1638.  
  1639.     aStart= fSelStart;
  1640.     if (!target || !*target) {  //search again
  1641.         if (fSearchRec.targ && (fSearchRec.step == -1 || fSearchRec.step == 1)) {  
  1642.             aSearchRec= fSearchRec;
  1643.             target= aSearchRec.targ;
  1644.             }        
  1645.         else {
  1646.             return kSearchNotFound;
  1647.             }
  1648.         }        
  1649.     else {
  1650.         limit= fLength; //StrLen(fBases);
  1651.         if (aStart < 0) aStart= 0;
  1652.         if (backwards) {
  1653.             if (!aStart) aStart= limit-1;
  1654.             else limit= aStart+1;  //? off-by-one ??
  1655.             }        
  1656.         else {
  1657.             if (aStart > limit) aStart= 0; 
  1658.             else limit -= aStart;
  1659.             }
  1660.         if (fSelBases > 0 && fSelBases < limit) limit= fSelBases;
  1661.         if (backwards) aSearchRec.step= -1;
  1662.         else aSearchRec.step= 1;
  1663.         aSearchRec.indx= aStart - aSearchRec.step;  //indx always in [0..handsize-1], except 1st is -1 
  1664.         aSearchRec.lim = limit + aSearchRec.step;  //lim always in [1..handsize] or 0=done, except 1st is hand+1 
  1665.         aSearchRec.targ= StrDup(target);
  1666.         }
  1667.  
  1668.     switch (fKind) {  
  1669.         case kRNA:
  1670.         case kDNA:
  1671.         case kNucleic: 
  1672.             NucleicSearch(fBases,target,aSearchRec); 
  1673.             break;
  1674.         default:
  1675.             ProteinSearch(fBases,target,aSearchRec);
  1676.             break;
  1677.         }
  1678.     
  1679.     fSearchRec= aSearchRec;
  1680.     if (aSearchRec.lim == 0) return kSearchNotFound;
  1681.     else return aSearchRec.indx;
  1682. }  
  1683.  
  1684.  
  1685.